/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

/*! \file Stochasticprocess.h
    \brief Stochasticprocess.h
*/

#include <src/miscellaneous/OS_basic.h>
#include <src/math/BoundaryCondition/BoundaryCondition.h>
#ifndef _OSStochastic_Processes_
#define _OSStochastic_Processes_


    //! prototype class for multi-dimensional stochastic process.
    /*! This class describes a stochastic process governed by
        \f[
        d\mathrm{x}_t = \mu(t, x_t)\mathrm{d}t
                      + \sigma(t, \mathrm{x}_t) \cdot d\mathrm{W}_t.
        \f]
    */
	template<class underlying_type, class matrix_type, OS_size Dim = 0, class time_type = OS_double>
	class OSStochasticProcess {
      public:
		typedef OSStochasticProcess											OS_Base;
		typedef underlying_type												underlying_type;
		typedef matrix_type													matrix_type;
		typedef time_type													time_type;
		typedef typename OS_type_helper<underlying_type>::type::base_value_type
																			base_value_type;
		enum :OS_size{Dim = Dim};
		enum :OS_size{Matrix_Dim = 2*Dim};

		OS_STATIC_CHECK( (2*Dim == Matrix_Dim) );

        //! discretization of a stochastic process over a given time interval
        virtual ~OSStochasticProcess() {}
        //! \name Stochastic process interface
        //@{
        //! returns the number of dimensions of the stochastic process
		OS_size size() const {return Dim +1;};
        //! returns the initial values of the state variables
        virtual underlying_type initial_conditions() const = 0;
        /*! \brief returns the rate of the equation, i.e.,
                   \f$ \mu(t, \mathrm{x}_t) \f$
        */
        virtual underlying_type mu(time_type t, const underlying_type& x) const = 0;
		template <class result_type, class argument_type>
		result_type mu_helper(time_type t, const argument_type & x) {
			result_type result(Dim);
			argument_type::const_iterator x_it = x.begin(), x_it_end=x.end();
			result_type::iterator result_it = result.begin();
			for (;x_it != x_it_end;++x_it,++result_it)
				*result_it = mu(t,*x_it);
			return result;
		};
        /*! \brief returns the diffusion part of the equation, i.e.
                   \f$ \sigma(t, \mathrm{x}_t) \f$
        */
        virtual matrix_type sigma(time_type t, const underlying_type& x) const = 0;
		template <class result_type, class argument_type>
		result_type sigma_helper(time_type t, const argument_type & x) {
			result_type result(x.size());
			argument_type::const_iterator x_it = x.begin(), x_it_end=x.end();
			result_type::iterator result_it = result.begin();
			for (;x_it != x_it_end;++x_it,++result_it)
				*result_it = sigma(t,*x_it);
			return result;
		};
		virtual OS_double compound(time_type t) const {
			OS_DYNAMIC_CHECK(false, "to implement in dervided class");
		};
      protected:
		  OSStochasticProcess(){};
    };

	template<class underlying_type>
	class drifted_brownian_process_base : public OSStochasticProcess<
		underlying_type,
		typename OS_Multi_Array<OS_Vector,2*OS_type_helper<underlying_type>::type::Dim,typename OS_type_helper<underlying_type>::type::base_value_type>::type
	>
	{

	public:
		typedef drifted_brownian_process_base OS_Base;

		drifted_brownian_process_base() 
			: rate_(Dim), volatility_(Dim), S0_(Dim)
		{};
		drifted_brownian_process_base(const underlying_type &rate, const underlying_type &volatility, const underlying_type &initial_conditions) 
			: rate_(rate), volatility_(volatility), S0_(initial_conditions)
		{
		};

		underlying_type initial_conditions() const {return S0_;};
		underlying_type volatility(time_type t, const underlying_type& x) const {return volatility_;};
		underlying_type rate(time_type t, const underlying_type& x) const {return rate_;};
		underlying_type mu(time_type t, const underlying_type& x) const {return rate_;};
		virtual OS_double compound(time_type t) const {
			OS_DYNAMIC_CHECK(false, "to implement in dervided class");
		};

		underlying_type		volatility_;
		underlying_type		rate_;
		underlying_type		S0_;
    };

    //! interface specialization for multi-dimensional drifted brownian motion.
    /*! This class describes a stochastic process governed by
        \f[
        d\mathrm{x}_t = \mu \mathrm{d}t
                      + \sigma \cdot d\mathrm{W}_t.
        \f]
    */

	template<class underlying_type, OS_size Dim = OS_type_helper<underlying_type>::type::Dim>
	class drifted_brownian_process : public drifted_brownian_process_base<underlying_type>
	{
      public:
		typedef drifted_brownian_process OS_Base;

		typedef drifted_brownian_process_base Base;
		drifted_brownian_process(OS_size size) : Base(size){};
		matrix_type sigma(time_type t, const underlying_type& x) const {return sigma_;};
		drifted_brownian_process(const underlying_type &rate, const underlying_type &volatility, const underlying_type &initial_conditions, const matrix_type& correlations) 
			: Base(rate,volatility,S0), sigma_(correlations)
		{
			OS_DYNAMIC_CHECK(false, "implement sigma_ = f(correlation,volatility)");
		};
		virtual OS_double compound(time_type t) const {
			OS_DYNAMIC_CHECK(false, "to implement");
		};
		matrix_type sigma_;
    };

    //! interface specialization for a one-dimensional drifted stochastic process.
    /*! This class describes a stochastic process governed by
        \f[
        d\mathrm{x}_t = \mu \mathrm{d}t
                      + \sigma d\mathrm{W}_t.
        \f]
		for a one dimensional random variable \f$\mathrm{x}_t\f$.
    */
	template<>
	class drifted_brownian_process<double,0> : public drifted_brownian_process_base<double>
	{
      public:
		typedef drifted_brownian_process_base Base;
		drifted_brownian_process(const underlying_type &rate, const underlying_type &volatility, const underlying_type &initial_conditions) 
			: Base(rate,volatility,initial_conditions) {};

		matrix_type sigma(time_type t, const underlying_type& x) const {return volatility_;};
		virtual OS_double compound(time_type t) const {
			OS_DYNAMIC_CHECK(false, "to implement");
		};
    };

    //! interface specialization for multi-dimensional drifted geometric motion.
    /*! This class describes a stochastic process governed by
        \f[
        d\mathrm{x}_t = \mu \mathrm{x}_t \mathrm{d}t
                      + \sigma \cdot \mathrm{x}_t d\mathrm{W}_t.
        \f]
    */
	template<class underlying_type = double, OS_size Dim = OS_type_helper<underlying_type>::type::Dim>
	class geometric_process : public drifted_brownian_process<underlying_type,Dim>
	{
      public:
		typedef drifted_brownian_process Base;
		typedef geometric_process OS_Base;
		geometric_process(OS_size size) : Base(size){};
		geometric_process(const underlying_type &rate, const underlying_type &volatility, const underlying_type &initial_conditions, const matrix_type& correlations) 
			: Base(rate,volatility,initial_conditions,correlations)
		{
		};
		underlying_type mu(time_type t, const underlying_type& x) const {underlying_type result(x); OS_Array_xor(result,rate_); return result;};
		matrix_type sigma(time_type t, const underlying_type& x) const { 
			return OS_Array_Tensor(sigma_,x);
		};
		underlying_type drift(time_type t, const underlying_type& x) const {return (rate_-.5*(volatility_^volatility_))*t;};
    };

    //! interface specialization for one-dimensional drifted geometric motion.
    /*! This class describes a stochastic process governed by
        \f[
        d\mathrm{x}_t = \mu \mathrm{x}_t \mathrm{d}t
                      + \sigma \mathrm{x}_t d\mathrm{W}_t.
        \f]
		for a one dimensional random variable \f$\mathrm{x}_t\f$.
    */
	template<>
	class geometric_process<double,0> : public drifted_brownian_process<double,0>
	{
      public:
	    typedef drifted_brownian_process Base;
		typedef geometric_process OS_Base;
		geometric_process(underlying_type rate, underlying_type volatility, underlying_type initial_conditions) : Base(rate,volatility,initial_conditions)
		  {};
		underlying_type mu(time_type t, const underlying_type& x) const {return rate_*x;};
		matrix_type sigma(time_type t, const underlying_type& x) const {return volatility_*x;};
		underlying_type drift(time_type t, const underlying_type& x) const {return rate_-.5*volatility_*volatility_;};
		virtual OS_double compound(time_type t) const {
			return std::exp(-rate_*t);
		};
    };

/* Cox Ingelsson Ross process Interface
*/
	template<class underlying_type = double>
	class CIR_process : public OSStochasticProcess<underlying_type,underlying_type>
	{
		OS_STATIC_CHECK(Dim == 0); // not known generalization in higher dimension
	public:
		typedef OSStochasticProcess Base;

		CIR_process() 
			: volatility_(0), kappa_(0), S0_(0), theta_(0)
		{};
		CIR_process(
			const underlying_type &sigma, 
			const underlying_type &kappa, 
			const underlying_type &theta,
			const underlying_type &S0)
			: sigma_(sigma), kappa_(kappa), theta_(theta),S0_(S0)
		{
			delta_ = 4.*kappa_*theta_/(sigma_*sigma_);
		};

		underlying_type initial_conditions() const {return S0_;};
		underlying_type sigma(time_type t, const underlying_type& x) const {return sigma_*std::sqrt(x);};
		underlying_type mu(time_type t, const underlying_type& x) const {return kappa_*(theta_-x);};
		virtual OS_double compound(time_type t) const {
            return std::exp(-B(0., t)*S0_);
		};
		underlying_type		sigma_;
		underlying_type		kappa_;
		underlying_type		theta_;
		underlying_type		S0_;
		underlying_type		delta_;
	protected: // taken from quantlib
    underlying_type A(time_type t, time_type T) const {
        underlying_type sigma2 = sigma_*sigma_;
        underlying_type h = std::sqrt(kappa_*kappa_ + 2.0*sigma2);
        underlying_type numerator = 2.0*h*std::exp(0.5*(kappa_+h)*(T-t));
        underlying_type denominator = 2.0*h + (kappa_+h)*(std::exp((T-t)*h) - 1.0);
        underlying_type value = std::log(numerator/denominator)*2.0*kappa_*theta_/sigma2;
        return std::exp(value);
    }

    underlying_type B(time_type t, time_type T) const {
        underlying_type h = std::sqrt(kappa_*kappa_ + 2.0*sigma_*sigma_);
        underlying_type temp = std::exp((T-t)*h) - 1.0;
        underlying_type numerator = 2.0*temp;
        underlying_type denominator = 2.0*h + (kappa_+h)*temp;
        underlying_type value = numerator/denominator;
        return value;
    }

    };
 
 
#endif
